1.1 LaTeX

R Markdown est encore plus puissant, on peut utilisier \(\LaTeX\), soit inline comme \(y_i=\beta_0 + x_i\beta_1\) ou comme equation: \[ \hat{\beta}=(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{y} \] vous pouvez même inclure des fichiers .tex entiers.

1.2 Figures

Imuages (ce n’est pas un plot) peuvent également être facilement incorporés avec ![](la_meilleure_image.png)

1.3 Knitting

Pour les documents de base, le yaml est ce qui est important:

rmdformats::robobook c’est seulement pour rendre le ficher légèrement plus joli

1.4 Pourquoi?

• Projets pour votre thèse de baccalauréat (Peut être entièrement reproductible)
• après votre baccalauréat
  • Industrie-1: Vous ferez beaucoup de rapports, souvent répétitifs, ce qui n’est pas très fun.. mais vous pouvez automatiser la plupart avec markdown
  • Industrie-2: Particulierement en tant que junior, vous aurez besoin de faire des recherches, Rmarkdown est parfait pour les présenter
  • Université: Si vous faites une maitrise - Vous devrez probablement coder et remettre votre code
  • Vous pouvez également utiliser python à partir des chunks de rmarkdown

2.1 SQL-Like

SQL (Structured Query Language) - votre meilleur ami si vous travaillez avec des données dans une entreprise, mais également utile à plus petite échelle (aka. pour vos projets).

L’exercice est basé sur l’article : “Where is the Land of Opportunity ? The Geography of Intergenerational Mobility in the United States” (Raj Chetty, Nathaniel Hendren, Patrick Kline et Emmanuel Saez ; Quarterly Journal of Economics 129(4) : 1553-1623, 2014).Il fait partie d’un plus grand projet disponible à l’adresse suivante : http://www.equality-of-opportunity.org/.

Nous commençons par charger quelques données échantillons. Notez que “kable” dans la dernière ligne rend seulement les graphiques jolis.

library(haven) # si le format est .dta
library(readxl) # se le format est .xls / .xlsx

mobility_data <- read_dta("data/mobility_across_cz.dta")
codebook <- read_excel("data/codebook.xls")

mobility_data %>% 
  head() %>% 
  kable()

Nous voyons qu’il y a des données manquantes et que les variables explicatives sont normalement numériques. Nous nous en occuperons plus tard. Pour une explication des variables, considérez le fichier codebook.xls.

select

En général, un jeu de données comporte de nombreuses colonnes. Mais parfois, nous n’avons besoin de considérer que certaines d’entre elles. La logique “select” vous permet de choisir un sous-ensemble des colonnes de données.

Il existe de nombreuses façons de le faire. Ci-dessous la démonstration pour sélectionner une seule colonne “e_rank_b” et une collection de colonnes c(“e_rank_b”, “taxrate”, “crime_violent”). C’est quoi la difference?

# Une seule colonne
facon_1 <- mobility_data$e_rank_b
facon_2 <- mobility_data[,'e_rank_b']
facon_3 <- mobility_data %>% select(e_rank_b)

# Plusieurs
facon_1_mul <- mobility_data[, c("e_rank_b", "taxrate", "crime_violent")]
facon_2_mul <- mobility_data %>% select(e_rank_b, taxrate, crime_violent)


facon_1 %>% head() 

## [1] 43.58841 39.64703 46.93892       NA 43.22648 53.63568

facon_1_mul %>% head()

## # A tibble: 6 × 3
##   e_rank_b taxrate crime_violent
##      <dbl>   <dbl>         <dbl>
## 1     43.6  0.0159        0.0483
## 2     39.6  0.0257        0.0316
## 3     46.9  0.0117        0.0512
## 4     NA    0.0217        0.0371
## 5     43.2  0.0452        0.0500
## 6     53.6  0.0267        0.0322

facon_2 %>% head() 

## # A tibble: 6 × 1
##   e_rank_b
##      <dbl>
## 1     43.6
## 2     39.6
## 3     46.9
## 4     NA  
## 5     43.2
## 6     53.6

facon_2_mul %>% head() 

## # A tibble: 6 × 3
##   e_rank_b taxrate crime_violent
##      <dbl>   <dbl>         <dbl>
## 1     43.6  0.0159        0.0483
## 2     39.6  0.0257        0.0316
## 3     46.9  0.0117        0.0512
## 4     NA    0.0217        0.0371
## 5     43.2  0.0452        0.0500
## 6     53.6  0.0267        0.0322

if-condition (filter)

Les commandes de “select” agissent sur les colonnes, les équivalents pour les lignes sont les “if-conditions” dans R, on appel ca les conditions “filter”. En reprenant l’exemple ci-dessus. Nous ne voulons que les observations où le e_rank_b est supérieur à 40. Quelle différence peut-on constater entre les deux méthodes ?

facon_1 <- mobility_data[mobility_data$e_rank_b > 40,
                         c("e_rank_b", "taxrate", "crime_violent")]

facon_2 <- mobility_data %>% 
  select(e_rank_b, taxrate, crime_violent) %>% 
  filter(e_rank_b > 40)

facon_1 %>% head()

## # A tibble: 6 × 3
##   e_rank_b taxrate crime_violent
##      <dbl>   <dbl>         <dbl>
## 1     43.6  0.0159        0.0483
## 2     46.9  0.0117        0.0512
## 3     NA   NA            NA     
## 4     43.2  0.0452        0.0500
## 5     53.6  0.0267        0.0322
## 6     45.6  0.0412        0.0267

facon_2 %>% head()

## # A tibble: 6 × 3
##   e_rank_b taxrate crime_violent
##      <dbl>   <dbl>         <dbl>
## 1     43.6  0.0159        0.0483
## 2     46.9  0.0117        0.0512
## 3     43.2  0.0452        0.0500
## 4     53.6  0.0267        0.0322
## 5     45.6  0.0412        0.0267
## 6     48.8  0.0222        0.0286

groupby

Parfois, nous sommes intéressés par des données agrégées. Avec la commande group_by, nous pouvons facilement créer des statistiques sommaires pour des sous-groupes de données.

Il existe des façons de faire exactement cela dans base R, mais ici nous nous concentrons sur la logique de dplyr, qui est plus intuitive.

# Calculer la moyenne de "e_rank_b" pour chaque ville

group_stat <- mobility_data %>% # Choisir jeu des données
  select(czname, e_rank_b) %>% # Choisir les colonnes
  group_by(czname) %>% # Specifier les groups
  summarise(
    mean_per_city = mean(e_rank_b, na.rm=T) # Specifier la statistique d'aggreagtion
  )

group_stat %>% head()

## # A tibble: 6 × 2
##   czname     mean_per_city
##   <chr>              <dbl>
## 1 Aberdeen            45.6
## 2 Aiken               38.5
## 3 Ainsworth           48.6
## 4 Albany              52.9
## 5 Alexandria          46.6
## 6 Allentown           61.3

(mutate)

Souvent, nous devons transformer les variables pour obtenir des informations significatives de notre analyse. Cela peut être fait facilement avec la logique “mutate”. L’utilisation de condtions if-else pour créer de nouvelles variables est extrêmement utile. Par exemple, nous pouvons diviser la mobilité sociale (e_rank_b) en deux catégories : faible et élevée.

data_tmp <- mobility_data
data_tmp <- data_tmp[!is.na(data_tmp$e_rank_b), ]
data_tmp$binary_mobility <- ifelse(data_tmp$e_rank_b > 45, 'High', 'Low')

ex_mutate_2 <- mobility_data %>% 
  filter(!is.na(e_rank_b)) %>% # PQ? 
  mutate(binary_mobility = if_else(e_rank_b > 45, 'High', 'Low')) %>% 
  select(binary_mobility, e_rank_b)


data_tmp[,c('binary_mobility', 'e_rank_b')] %>% head()

## # A tibble: 6 × 2
##   binary_mobility e_rank_b
##   <chr>              <dbl>
## 1 Low                 43.6
## 2 Low                 39.6
## 3 High                46.9
## 4 Low                 43.2
## 5 High                53.6
## 6 Low                 38.1

ex_mutate_2 %>% head()

## # A tibble: 6 × 2
##   binary_mobility e_rank_b
##   <chr>              <dbl>
## 1 Low                 43.6
## 2 Low                 39.6
## 3 High                46.9
## 4 Low                 43.2
## 5 High                53.6
## 6 Low                 38.1

2.2 Exercises

1. Transformer la variable “e_rank_b” en “inférieure à 50” et “supérieure à 50”. Calculez la moyenne de “crime_violent” pour chaque groupe. Quelles sont vos conclusions?
2. Transformer la variable “frac_middleclass” en “inférieure à 0.5” et “supérieure à 0.5”. Pareil pour la variable “eitc_exposure”. Calculer la moyenne du taux d’abandon scolaire au lycée “dropout_r” pour chaque combinaison des deux nouvelles variables. Conseil : utilisez group_by
3. Proposer une façon de traiter les valeurs “NA”. Pourquoi pensez-vous que c’est une bonne idée ?

3.1 Plots

Les plots/figures peuvent facilement être intégrées dans markdown. Tout comme une sortie de tableau, on peut les structurer légèrement. Par exemple, nous pouvons ajouter une légende ou ajouter plusieurs graphiques côte à côte.

Par exemple, commençons par un histogramme:

# Histogram simple
hist(mobility_data$e_rank_b)

?hist

?hist

plot(density(data_tmp$e_rank_b)) # Pourquoi changer les données?
plot(density(data_tmp$e_rank_b), col='red')

Il existe aussi un package sympa de l’écosystème tidyverse (ggplot2)

plot_1 <- mobility_data %>% 
  ggplot() + 
  geom_histogram(aes(x=e_rank_b))

plot_2 <- mobility_data %>% 
  ggplot() + 
  geom_histogram(aes(x=crime_violent))

grid.arrange(plot_1, plot_2, ncol=2)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 19 rows containing non-finite values (stat_bin).

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 17 rows containing non-finite values (stat_bin).

Nous pouvons également construire des figures avec plusieurs informations. Par exemple, un histogramme avec une ligne de densité. (Attention toutefois aux axes!)

# Base R
hist(data_tmp$e_rank_b, probability = T)
lines(density(data_tmp$e_rank_b), col='red')

# GGplot
data_tmp %>% 
  ggplot() + 
  geom_histogram(aes(x=e_rank_b, y=..density..)) + 
  geom_density(aes(x=e_rank_b), color='red')

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

3.2 Exercises

Conseil : si vous n’arrivez pas à trouver la bonne commande pour une figure, les solutions, utilisez google!

1. Montrer la relation entre le taux d’obtention d’un diplôme universitaire et la mobilité socialeà l’aide d’un diagramme de dispersion (scatterplot) et d’une ligne de régression. Conseil: vous pouvez trouver une solution si vous googlez “r plot with regression line”

2. Existe-t-il une relation entre les dépenses des gouvernements locaux par habitant (subcty_total_expenditure_pc) et le taux d’imposition (taxrate)? Pouvez-vous utiliser les données directement ? Conseil : pour mettre les données à l’échelle, vous pouvez toujours essayer d’utiliser log, sqrt, etc.

3. Montrez graphiquement qu’il existe une ville où la criminalité violente est exceptionnellement élevée (‘crime_violent’). Que pouvez-vous faire avec de telles données (outlier ou pas)?

4.1 Tables

Les tableaux peuvent facilement être incorporés dans R. Vous pouvez soit les imprimer directement, soit utiliser kable (du paquet knitr - plus jolie). Bien sûr, nous pouvons également combiner ce que nous avons vu dans les sections précédentes (en particulier la manipulation des données)

# Impression de base
mobility_data %>% head()

## # A tibble: 6 × 39
##   index    cz czname         e_rank_b cs_race_bla cs_race_theil_20… cs00_seg_inc
##   <dbl> <dbl> <chr>             <dbl>       <dbl>             <dbl>        <dbl>
## 1     0   100 Johnson City       43.6    0.0694             0.0577        0.0409
## 2     1   200 Morristown         39.6    0.445              0.0750       -0.0259
## 3     2   301 Middlesborough     46.9    0.0151             0.110         0.0954
## 4     3   302 Knoxville          NA     -0.000411          -0.00446       0.0105
## 5     4   401 Winston-Salem      43.2   -0.0193             0.263         0.0289
## 6     5   402 Martinsville       53.6    0.00563            0.00699       0.0221
## # … with 32 more variables: cs00_seg_inc_pov25 <dbl>, cs00_seg_inc_aff75 <dbl>,
## #   frac_traveltime_lt15 <dbl>, hhinc00 <dbl>, gini <dbl>,
## #   inc_share_1perc <dbl>, gini99 <dbl>, frac_middleclass <dbl>, taxrate <dbl>,
## #   subcty_total_expenditure_pc <dbl>, tax_st_diff_top20 <dbl>,
## #   eitc_exposure <dbl>, ccd_exp_tot <dbl>, ccd_pup_tch_ratio <dbl>,
## #   score_r <dbl>, dropout_r <dbl>, num_inst_pc <dbl>, tuition <dbl>,
## #   gradrate_r <dbl>, cs_labforce <dbl>, cs_elf_ind_man <dbl>, …

# Avec Kable
mobility_data %>% head() %>% kable() # c'est comme kable(head(mobilitiy_data))

ici nous pouvons manipuler légèrement nos données

# Example avec manipulation des données 
# Décrire ce qui se passe ici
mobility_data %>% 
  filter(!is.na(e_rank_b)) %>% 
  mutate(decile = ntile(e_rank_b, 10)) %>% 
  select(decile, e_rank_b, mig_outflow) %>% 
  group_by(decile) %>% 
  summarise(average_outflow = mean(mig_outflow)) %>% 
  kable()

4.2 Distributions en R

L’un des exercices les plus importants lors de la modélisation statistique consiste à vérifier les hypothèses du modèle. La logique est toujours la même.

q-distribution (quantiles, ex: qnorm, qpois, qgamma etc.) Quel est le 99ème quantile d’une gamma avec shape=2 et rate=1?

qgamma(0.99, shape=2,rate=1)

## [1] 6.638352

d-distribution (densité, ex: dnorm, dpois) Tracez visuellement le 99ème quantile d’une gamma avec shape=2 et rate=1?

gamma_dist <- dgamma(seq(0,8,0.1), shape=2,rate=1)

plot(x=seq(0,8,0.1), 
     y=gamma_dist,
     type='l')
abline(v=qgamma(0.99, shape=2,rate=1), 
       col='red', lty=2)

p-distribution (densité cumulative, ex: pnorm, ppois)

Supposons que \(X \sim Gamma(2,1)\). Quel pourcentage de la masse totale se trouve à droite de 5 ?

1 - pgamma(5, shape=2, rate=1)

## [1] 0.04042768

r-distribution (échantillon, ex: rnorm, rpois)

Simulez 100 observations à partir d’un gamma avec shape=2, rate=1 et comparez-les à la distribution théorique.

set.seed(42)
sims_gamma <- rgamma(100, shape=2, rate=1)

plot(density(sims_gamma), ylim=c(0,0.38), 
     xlim=c(-1, 7), 
     main='Comparaison empirique vs. theoretique')
lines(x=seq(0,8,0.1), 
      y=gamma_dist, col='red')

# A la main
# Sans transformation
inc <- (as.numeric(mobility_data$hhinc00))
inc <- (inc - mean(inc)) / sd(inc)

quantiles_norm <- qnorm(seq(0,1,0.001))
quantiles_empirical <- quantile(inc, seq(0,1,0.001))

plot(quantiles_norm,
     quantiles_empirical, 
     main='Normal Q-Q Plot', 
     ylab='Sample Quantiles',  
     xlab='Theoretical Quantiles')
abline(a = 0, b = 1, col = "red")

# Avec transformation
# A la main
# Sans transformation
log_inc <- log(as.numeric(mobility_data$hhinc00))
log_inc <- (log_inc - mean(log_inc)) / sd(log_inc)

quantiles_norm <- qnorm(seq(0,1,0.001))
quantiles_empirical <- quantile(log_inc, seq(0,1,0.001))

plot(quantiles_norm,
     quantiles_empirical, 
     main='Normal Q-Q Plot (Log transformed)', 
     ylab='Sample Quantiles',  
     xlab='Theoretical Quantiles')
abline(a = 0, b = 1, col = "red")

qqnorm(log(as.numeric(mobility_data$hhinc00)))
qqline(log(as.numeric(mobility_data$hhinc00)),
       col = 2,lwd=2,lty=2)

4.3 Exercises

1. Supposons que \(X \sim \mathcal{N}(0,10)\). Echantillonnez 10 observations et calculez la moyenne, que nous appelons Y. Répétez l’opération ci-dessus 1000 fois et calculez la moyenne et l’écart-type de Y. A quoi vous attendez-vous et pourquoi, qu’observez-vous ?

2. Créez un tableau ou qui montre le nombre d’observations manquantes pour chaque variable (colonne dans l’ensemble de données).

3. Analysez les variables tuition et num_inst_pc et discutez de ce que vous voyez. Quelles sont vos conclusions ? S’agit-il d’un problème?

4. Créez un tableau montrant le coefficient de corrélation entre la mobilité sociale e_rank_b et toutes les autres variables. Évaluez vos résultats.

5. Séparez la variable en “faible”, “moyen” et “élevé”. Justifiez votre choix. Comparez chaque niveau avec les variables gini et cs_born_foreign. Interprétez et créez un tableau.

5.1 Charger et preparer les données

# Charger l'essentiel
library(tidyverse)
library(knitr)
library(gridExtra)
library(haven) 
library(readxl)

mobility_data <- read_dta("data/mobility_across_cz.dta")
codebook <- read_excel("data/codebook.xls")

mobility_data %>% 
  head() %>% 
  kable()

str(mobility_data)

## tibble [495 × 39] (S3: tbl_df/tbl/data.frame)
##  $ index                      : num [1:495] 0 1 2 3 4 5 6 7 8 9 ...
##   ..- attr(*, "format.stata")= chr "%12.0g"
##  $ cz                         : num [1:495] 100 200 301 302 401 402 500 601 602 700 ...
##   ..- attr(*, "format.stata")= chr "%12.0g"
##  $ czname                     : chr [1:495] "Johnson City" "Morristown" "Middlesborough" "Knoxville" ...
##   ..- attr(*, "format.stata")= chr "%22s"
##  $ e_rank_b                   : num [1:495] 43.6 39.6 46.9 NA 43.2 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_race_bla                : num [1:495] 0.069393 0.445318 0.015122 -0.000411 -0.019341 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_race_theil_2000         : num [1:495] 0.05769 0.07497 0.11006 -0.00446 0.26349 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs00_seg_inc               : num [1:495] 0.0409 -0.0259 0.0954 0.0105 0.0289 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs00_seg_inc_pov25         : num [1:495] -0.0171 0.0353 0.0139 0.0545 0.0478 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs00_seg_inc_aff75         : num [1:495] 0.02664 0.0089 0.05273 0.00763 0.06856 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ frac_traveltime_lt15       : num [1:495] 0.384 0.354 0.479 0.661 0.406 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ hhinc00                    : num [1:495] 38210 25123 32701 36524 30403 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ gini                       : num [1:495] 0.419 0.487 0.307 0.241 0.451 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ inc_share_1perc            : num [1:495] 5.7 11.82 11.81 NA 9.97 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ gini99                     : num [1:495] 0.253 0.393 0.208 NA 0.348 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ frac_middleclass           : num [1:495] 0.584 0.356 0.668 NA 0.538 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ taxrate                    : num [1:495] 0.0159 0.0257 0.0117 0.0217 0.0452 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ subcty_total_expenditure_pc: num [1:495] 1405 1569 3116 2514 2138 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ tax_st_diff_top20          : num [1:495] 1.07949 0.01987 -0.00462 0.20946 0.0028 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ eitc_exposure              : num [1:495] 0.000078 -0.012866 0.034301 1.929972 -0.003527 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ ccd_exp_tot                : num [1:495] 6.58 4.41 6.61 7.89 7.05 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ ccd_pup_tch_ratio          : num [1:495] 18 17.8 17.9 12.6 20.9 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ score_r                    : num [1:495] NA -16.02 -9.14 NA -10.45 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ dropout_r                  : num [1:495] NA 0.0562 0.0246 NA NA ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ num_inst_pc                : num [1:495] 0.0554 0.0379 0.0635 NA 0.0419 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ tuition                    : num [1:495] 1032 970 1395 NA 2524 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ gradrate_r                 : num [1:495] -0.221 -0.0416 -0.1717 NA 0.1363 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_labforce                : num [1:495] 0.653 0.527 0.707 0.623 0.636 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_elf_ind_man             : num [1:495] 0.06 0.2087 0.0363 0.0685 0.1408 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ d_tradeusch_pw_1990        : num [1:495] NA 0.767 0.0275 0.0173 0.8333 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ frac_worked1416            : num [1:495] 0.0149 0.0224 0.0275 NA 0.0239 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ mig_inflow                 : num [1:495] 0.0417 0.0182 0.031 NA 0.058 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ mig_outflow                : num [1:495] 0.0559 0.0138 0.0422 NA 0.0682 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_born_foreign            : num [1:495] 0.1605 0.0232 0.1318 0.0253 0.1903 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ scap_ski90pcm              : num [1:495] NA -0.613 -2.064 2.343 -0.694 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ rel_tot                    : num [1:495] 0.536 0.647 0.422 0.688 0.387 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ crime_violent              : num [1:495] 0.0483 0.0316 0.0512 0.0371 0.05 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_fam_wkidsinglemom       : num [1:495] 0.242 0.35 0.183 0.129 0.222 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_divorced                : num [1:495] 0.139 0.104 0.14 0.113 0.115 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"
##  $ cs_married                 : num [1:495] 0.549 0.543 0.622 0.71 0.57 ...
##   ..- attr(*, "format.stata")= chr "%10.0g"

Le format est celui de stata, mais la plupart des variables sont numériques, cela semble correct.

# Verifier données manquantes
for (column_ in names(mobility_data)){
  print(
    paste(
      column_, 
      sum(is.na(mobility_data[, column_])), 
      sep=': '
    )
  )
}

## [1] "index: 0"
## [1] "cz: 0"
## [1] "czname: 0"
## [1] "e_rank_b: 19"
## [1] "cs_race_bla: 0"
## [1] "cs_race_theil_2000: 0"
## [1] "cs00_seg_inc: 0"
## [1] "cs00_seg_inc_pov25: 0"
## [1] "cs00_seg_inc_aff75: 0"
## [1] "frac_traveltime_lt15: 0"
## [1] "hhinc00: 0"
## [1] "gini: 0"
## [1] "inc_share_1perc: 19"
## [1] "gini99: 19"
## [1] "frac_middleclass: 19"
## [1] "taxrate: 1"
## [1] "subcty_total_expenditure_pc: 2"
## [1] "tax_st_diff_top20: 0"
## [1] "eitc_exposure: 0"
## [1] "ccd_exp_tot: 7"
## [1] "ccd_pup_tch_ratio: 20"
## [1] "score_r: 21"
## [1] "dropout_r: 89"
## [1] "num_inst_pc: 100"
## [1] "tuition: 104"
## [1] "gradrate_r: 107"
## [1] "cs_labforce: 0"
## [1] "cs_elf_ind_man: 0"
## [1] "d_tradeusch_pw_1990: 14"
## [1] "frac_worked1416: 19"
## [1] "mig_inflow: 10"
## [1] "mig_outflow: 10"
## [1] "cs_born_foreign: 0"
## [1] "scap_ski90pcm: 14"
## [1] "rel_tot: 0"
## [1] "crime_violent: 17"
## [1] "cs_fam_wkidsinglemom: 0"
## [1] "cs_divorced: 0"
## [1] "cs_married: 0"

Le nombre de points de données manquants ne semble pas non plus trop inquiétant. Il pourrait y avoir une certaine corrélation entre les points manquants mais nous pouvons ignorer cela pour le moment.

5.2 Renommer une variable

# renommer avec espace
print(names(mobility_data))

##  [1] "index"                       "cz"                         
##  [3] "czname"                      "e_rank_b"                   
##  [5] "cs_race_bla"                 "cs_race_theil_2000"         
##  [7] "cs00_seg_inc"                "cs00_seg_inc_pov25"         
##  [9] "cs00_seg_inc_aff75"          "frac_traveltime_lt15"       
## [11] "hhinc00"                     "gini"                       
## [13] "inc_share_1perc"             "gini99"                     
## [15] "frac_middleclass"            "taxrate"                    
## [17] "subcty_total_expenditure_pc" "tax_st_diff_top20"          
## [19] "eitc_exposure"               "ccd_exp_tot"                
## [21] "ccd_pup_tch_ratio"           "score_r"                    
## [23] "dropout_r"                   "num_inst_pc"                
## [25] "tuition"                     "gradrate_r"                 
## [27] "cs_labforce"                 "cs_elf_ind_man"             
## [29] "d_tradeusch_pw_1990"         "frac_worked1416"            
## [31] "mig_inflow"                  "mig_outflow"                
## [33] "cs_born_foreign"             "scap_ski90pcm"              
## [35] "rel_tot"                     "crime_violent"              
## [37] "cs_fam_wkidsinglemom"        "cs_divorced"                
## [39] "cs_married"

names(mobility_data)[4] <- "mobilite intergen"

# mobility_data$`mobilite intergen` -> On a toujours besoin des symboles `. Cela semble encombrant. 
# Essayez d'utiliser des tirets bas (_) ou des points (.)

names(mobility_data)[4] <- "mobilite_intergen"
# mobility_data$mobilite_intergen -> C'est mieux

5.3 Analyse visuelle de la mobilité

Nous allons analyser la densité en utilisant le lissage par noyau, un histogramme et la fonction de distribution cumulative.

hist(mobility_data$mobilite_intergen, 
     main='Intergenerational Mobility Scores', 
     xlab='Value', 
     ylab='Density', 
     probability = T, 
     breaks = seq(27, 67, 3))

lines(density(mobility_data$mobilite_intergen[!is.na(mobility_data$mobilite_intergen)]), 
      col='red', 
      lwd=2, 
      lty=2)

legend(55,0.07,
       legend=c("Density"),
       col=c("red"),
       lty=2,
       cex=0.8)

quantiles_ <- quantile(mobility_data$mobilite_intergen, 
                       probs = seq(0,1,0.01), 
                       na.rm = T)

plot(x=quantiles_, 
     y=seq(0,1,0.01), 
     type='l', 
     lwd=2, 
     main='ECDF', 
     ylab='', 
     xlab='Value for Mobility')

5.4 Analyse visuelle: Manufacturing vs Mobilité

plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

# Model simple pour la relation lineaire
summary(lm(mobilite_intergen ~ cs_elf_ind_man, data=mobility_data))

## 
## Call:
## lm(formula = mobilite_intergen ~ cs_elf_ind_man, data = mobility_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5226  -3.6620  -0.8982   3.1683  17.9719 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     46.9024     0.5867  79.946  < 2e-16 ***
## cs_elf_ind_man -17.7064     3.1306  -5.656 2.68e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.482 on 474 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.06322,    Adjusted R-squared:  0.06125 
## F-statistic: 31.99 on 1 and 474 DF,  p-value: 2.683e-08

plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

abline(a=46.9024,
       b=-17.7064,
       col='red', 
       lwd=2, 
       lty=2)

legend(0.35,60,
       legend=c("Linear Fit"),
       col=c("red"),
       lty=2,
       cex=0.8)

5.5 Analyse descriptif (secteur manufacturier)

Nous pouvons trouver de bons seuils en regardant le graphique ci-dessus ou en utilisant les quantiles. Voici l’approche “visuelle”. Lorsque nous coupons une variable, nous pouvons introduire des non-linéarités dans la relation (voir l’exercice suivant)

mobility_data$manuf_cut <- cut(mobility_data$cs_elf_ind_man, 
                               breaks = c(0,0.1,0.3, 1), 
                               labels = c('low', 'mid', 'high'))

mobility_data %>% 
  group_by(manuf_cut) %>% 
  summarise(
    mean_labforce = mean(cs_labforce, na.rm=T), 
    mean_growth_imports = mean(d_tradeusch_pw_1990, na.rm=T), 
    mean_mig_inflow = mean(mig_outflow, na.rm=T)
  )

## # A tibble: 3 × 4
##   manuf_cut mean_labforce mean_growth_imports mean_mig_inflow
##   <fct>             <dbl>               <dbl>           <dbl>
## 1 low               0.633               0.289          0.0458
## 2 mid               0.641               1.30           0.0427
## 3 high              0.667               2.73           0.0399

5.6 Model simple

model_lineaire_simple = lm(mobilite_intergen ~ manuf_cut, 
                           data=mobility_data)

summary(model_lineaire_simple)

## 
## Call:
## lm(formula = mobilite_intergen ~ manuf_cut, data = mobility_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.271  -3.847  -0.658   2.952  17.605 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    45.9575     0.5388  85.293  < 2e-16 ***
## manuf_cutmid   -2.3677     0.6141  -3.856 0.000131 ***
## manuf_cuthigh  -6.3979     1.2444  -5.142    4e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.495 on 473 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.06088,    Adjusted R-squared:  0.05691 
## F-statistic: 15.33 on 2 and 473 DF,  p-value: 3.534e-07

plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

x <- seq(0, 0.45, 0.01)
fx <- 45.9575 + (x > 0 & x <=0.1) * 0 +
      (x > 0.1 & x <= 0.3) * -2.3677 +
      (x > 0.3 & x <= 2.119) *  -6.3979
lines(x, fx, 
      col='red', 
      lty=2, 
      lwd=2)

plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

predictions_simple <- predict(model_lineaire_simple, 
                              mobility_data)

points(mobility_data$cs_elf_ind_man, 
       predictions_simple, 
       col='red')

model_lineaire_interaction = lm(mobilite_intergen ~ manuf_cut + cs_elf_ind_man + manuf_cut*cs_elf_ind_man, 
                           data=mobility_data)

summary(model_lineaire_interaction)

## 
## Call:
## lm(formula = mobilite_intergen ~ manuf_cut + cs_elf_ind_man + 
##     manuf_cut * cs_elf_ind_man, data = mobility_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.9773  -3.6249  -0.6577   2.7890  17.3450 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    42.921      1.983  21.640   <2e-16 ***
## manuf_cutmid                    3.175      2.252   1.410   0.1592    
## manuf_cuthigh                  -5.908     10.728  -0.551   0.5821    
## cs_elf_ind_man                 44.485     27.973   1.590   0.1124    
## manuf_cutmid:cs_elf_ind_man   -57.912     28.507  -2.031   0.0428 *  
## manuf_cuthigh:cs_elf_ind_man  -37.336     40.598  -0.920   0.3582    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.463 on 470 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.07768,    Adjusted R-squared:  0.06787 
## F-statistic: 7.917 on 5 and 470 DF,  p-value: 3.562e-07

plot(
  mobility_data$cs_elf_ind_man,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

predictions_interact <- predict(model_lineaire_interaction, 
                              mobility_data)

points(mobility_data$cs_elf_ind_man, 
       predictions_interact, 
       col='red')

5.7 Model “Classe moyenne”

plot(x=mobility_data$frac_middleclass, 
     y=mobility_data$mobilite_intergen, 
     main='Middle Class and Intergen. Mobility', 
     ylab='Intergen. Mobility', 
     xlab='Fraction Middle class')

model_middle_class <- lm(mobilite_intergen ~ frac_middleclass, 
                         data=mobility_data)

summary(model_middle_class)

## 
## Call:
## lm(formula = mobilite_intergen ~ frac_middleclass, data = mobility_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.9940  -2.9766  -0.5619   2.7327  15.9744 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        18.986      1.381   13.74   <2e-16 ***
## frac_middleclass   45.411      2.491   18.23   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.343 on 474 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.4121, Adjusted R-squared:  0.4109 
## F-statistic: 332.3 on 1 and 474 DF,  p-value: < 2.2e-16

plot(
  mobility_data$frac_middleclass,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

predictions_middleclass <- predict(model_middle_class, 
                              mobility_data)

points(mobility_data$frac_middleclass,
       predictions_middleclass, 
       col='red')

# Analyse des residues

residuals_ <- mobility_data$mobilite_intergen - predictions_middleclass
plot(mobility_data$frac_middleclass, 
     residuals_, 
     ylab='Residuals', 
     xlab='Middle class')

mobility_data_nomissings <- mobility_data[!is.na(mobility_data$frac_middleclass), ]

model_middle_class_poly <- lm(mobilite_intergen ~ poly(frac_middleclass, 2), 
                              data=mobility_data_nomissings)

summary(model_middle_class_poly)

## 
## Call:
## lm(formula = mobilite_intergen ~ poly(frac_middleclass, 2), data = mobility_data_nomissings)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.0057  -2.9207  -0.3494   2.4430  15.5096 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 43.9039     0.1881 233.464  < 2e-16 ***
## poly(frac_middleclass, 2)1  79.1684     4.1029  19.296  < 2e-16 ***
## poly(frac_middleclass, 2)2  31.2715     4.1029   7.622 1.38e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.103 on 473 degrees of freedom
## Multiple R-squared:  0.4764, Adjusted R-squared:  0.4742 
## F-statistic: 215.2 on 2 and 473 DF,  p-value: < 2.2e-16

plot(
  mobility_data$frac_middleclass,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

predictions_middleclass_poly <- predict(model_middle_class_poly,
                                        mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ])

lines(mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ]$frac_middleclass,
       predictions_middleclass_poly, 
       col='red')

5.8 Challenge

# List pour les resultats
results_list <- list()

# Loop
for (length_ in seq(1,10,1)){
  reg_ <- lm(mobilite_intergen ~ poly(frac_middleclass, length_), 
             data=mobility_data_nomissings)
  
  print(paste('R-Squared: ', summary(reg_)$r.squared))
  
  predictions_ <- predict(reg_,
                          mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ])
  
  results_list[[length_]] <- predictions_

}

## [1] "R-Squared:  0.412134418340672"
## [1] "R-Squared:  0.476437415437812"
## [1] "R-Squared:  0.477700187650047"
## [1] "R-Squared:  0.477727823275877"
## [1] "R-Squared:  0.478923076244006"
## [1] "R-Squared:  0.48017980394324"
## [1] "R-Squared:  0.48426795442976"
## [1] "R-Squared:  0.487546350228303"
## [1] "R-Squared:  0.487606950061351"
## [1] "R-Squared:  0.489193655479349"

plot(
  mobility_data$frac_middleclass,
  mobility_data$mobilite_intergen, 
  main='% Manufacturing vs Mobility', 
  ylab='Intergen. Mob', 
  xlab='Fraction Working in Manufacturing'
) 

for(i in seq(1,10,1)){
  
  predictions_middleclass_poly <- results_list[[i]]
  
  lines(mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ]$frac_middleclass,
         predictions_middleclass_poly, 
         col='red')
  
}